#
#The Following R-Code estimates OP, ZiOP, and ZiOPC models of political violence,
#using datasets similar to those used by Besley and Persson (2009)
#
#These results correspond to the models reported in Table 1 of our main paper
#
#
#All OP, ZiOP, and ZiOPC models are then compared with AIC statistics, 
#Vuong tests, and likelihood ratio tests.
#
#
#

#clear memory
rm( list=ls() ) 						      

#install.packages
library(mvtnorm)
library(car)                                                         
library(scatterplot3d)
library(plotrix)
library(Hmisc)
library(Zelig)
library(foreign, pos=4)
library(gplots)

#set the seed
set.seed(2)

#read in formatted replication dataset
BP.data <- read.dta("bp_exact_for_analysis.dta", 
  convert.dates=TRUE, convert.factors=TRUE, missing.type=TRUE, 
  convert.underscore=TRUE, warn.missing.labels=TRUE)

#summarize dataset
summary(BP.data)

#omit observations with missing data from dataset
BP.data<-na.omit(BP.data)

#summarize new dataset to ensure missignness has been removed
summary(BP.data)


#Begin Analysis#

############################################################################
##############################OP Model######################################
############################################################################

#This first model estimates an OP model of political violence comparable to B&P (2009)

#write a function to estimate an OP model 
OP <- function(b,data) {					      
  y<-data[,8] 			#political violence                                                        
  x1<-data[,4]			#lgdppc                                           
  x2<-data[,9]			#parliament
  x3<-data[,7]			#disaster
  x4<-data[,5]			#oil
  x5<-data[,6]			#primary                                               	                                             					     	      
  z<-cbind(x1,x2,x3,x4,x5)	#covariates
  tau<- rep(0,4) 		                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- b[1]+exp(b[2])
  tau[4]<- (Inf)
  n=nrow(data)							      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      
  for(i in 1:n){					              					      				      
        	G<-b[3:7]
       		ZG  <- G %*% z[i,]		              
                if((y[i]==0)){llik[i]<-log((pnorm(tau[2]-ZG)))}
                else if((y[i]==1)){llik[i]<-log((pnorm(tau[3]-ZG) - pnorm(tau[2]-ZG)))}
                else if((y[i]==2)) {llik[i]<-log((1 - pnorm(tau[3]-ZG)))}     
		}
        llik<--1*sum(llik)					      
    	return(llik)
}								      

#give the OP model resonable starting parameters
b<-rbind(-1, 0.3, -0.2, -0.5, 0.2,.9,-.4)

#optimize the OP model
output.OP<-nlm(f=OP, p=b, data=BP.data, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.OP)

#recover the variance covariance matrix
vcv<-solve(output.OP$hessian)

#show the variance covariance matrix to recover the standard errors
vcv 





############################################################################
##########################Main ZIOP Model###################################
############################################################################

#This model estimates our main ZiOP model of political violence

#write a function to estimate a ZiOP model 
ZIOP <- function(b,data) {					      
  y<-data[,8]                                                         
  x1<-data[,4]			#lgdppc                                           
  x2<-data[,9]			#parliament
  x3<-data[,7]			#disaster
  x4<-data[,5]			#oil
  x5<-data[,6]			#primary
  x<-cbind(1,x1,x2)		#inflation stage covariates 			     	      
  z<-cbind(x1,x2,x3,x4,x5)	#outcome stage covariates 

  tau<- rep(0,4) 			                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- b[1]+exp(b[2])
  tau[4]<- (Inf)
  n=nrow(data)							      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      
  for(i in 1:n){					              					      
        	B<-b[3:5]						      
 		XB  <- B %*% x[i,] 					      
        	G<-b[6:10]
       		ZG  <- G %*% z[i,]		              
                if((y[i]==0)){llik[i]<-log( (1-pnorm(XB))+ (pnorm(XB)) * (pnorm(tau[2]-ZG)) )}
                else if((y[i]==1)){llik[i]<-log( (pnorm(XB)) * (pnorm(tau[3]-ZG) - pnorm(tau[2]-ZG)))}
                else if((y[i]==2)) {llik[i]<-log( (pnorm(XB)) * (1-pnorm(tau[3]-ZG)))}     
		}
        llik<--1*sum(llik)					      
    	return(llik)
}


#give the ZiOP model resonable starting parameters
b<-rbind( -1.31, .32, 2.5, -.21,.2, -0.2, -0.4, 0.2,.9,-.4)

#optimize the ZiOP model
output.ZIOP<-nlm(f=ZIOP, p=b, iterlim=500, data=BP.data, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.ZIOP)

#recover the variance covariance matrix
vcv<-solve(output.ZIOP$hessian)

#show the variance covariance matrix to recover the standard errors
vcv  






############################################################################
##########################Small ZIOPC Model#################################
############################################################################

#This model estimates our smaller ZiOPC model of political violence

#write a function to estimate a ZiOPC model 
ZIOPC <- function(b,data) {					      
  y<-data[,8]			#political violence                                                         
  x1<-data[,4]			#lgdppc                                           
  x2<-data[,9]			#parliament
  x3<-data[,7]			#disaster
  x4<-data[,5]			#oil
  x5<-data[,6]			#primary
  x<-cbind(1,x1,x2)		#inflation stage covariates 			     	      
  z<-cbind(x1,x2,x3,x4,x5)	#outcome stage covariates 


  tau<- rep(0,4) 			                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- exp(b[2])+b[1]
  tau[4]<- (Inf)
  n=nrow(data)								      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      
  for(i in 1:n){					              					      
        	B<-b[3:5]						      
 		XB  <- B %*% x[i,] 				      
        	G<-b[6:10]
       		ZG  <- G %*% z[i,]   
                rho<-b[11] 
                sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
                negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)		              
                if((y[i]==0)){llik[i]<-log((1-pnorm(XB))+ pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
                else if((y[i]==1)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[3]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
                else if((y[i]==2)) {llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(ZG-tau[3])),mean=c(0,0),sigma=sigma))}     
		}
        llik<--1*sum(llik)					      
    	return(llik)
}	

############################################################################

#give the ZiOPC model resonable starting parameters
b<-rbind(  0.77,  0.90, 18.78, -2,.2, 0.04, -0.09,  0.26,  1.70, -0.42,-.1)

#optimize the ZiOPC model
output.ZIOPC<-nlm(f=ZIOPC, p=b, iterlim=500, data=BP.data, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.ZIOPC)

#recover the variance covariance matrix
vcv<-solve(output.ZIOPC$hessian)

#show the variance covariance matrix to recover the standard errors
vcv 





############################################################################
##########################Large ZIOPC Model#################################
############################################################################

#This model estimates our larger ZiOPC model of political violence

#write a function to estimate a ZiOPC model 
ZIOPC <- function(b,data) {					      
  y<-data[,8]			#political violence                                                         
  x1<-data[,4]			#lgdppc                                           
  x2<-data[,9]			#parliament
  x3<-data[,7]			#disaster
  x4<-data[,5]			#oil
  x5<-data[,6]			#primary
  x<-cbind(1,x1,x2,x3,x4,x5)	#inflation stage covariates 			     	      
  z<-cbind(x1,x2,x3,x4,x5)	#outcome stage covariates 

  tau<- rep(0,4) 			                              
  tau[1]<--(Inf)
  tau[2]<- b[1]
  tau[3]<- exp(b[2])+b[1]
  tau[4]<- (Inf)
  n=nrow(data)								      							      
  llik <- matrix(0, nrow=n, ncol = 1)   			      					      
  for(i in 1:n){					              					      
        	B<-b[3:8]						      
 		XB  <- B %*% x[i,] 				      
        	G<-b[9:13]
       		ZG  <- G %*% z[i,]   
                rho<-b[14] 
                sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
                negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)		              
                if((y[i]==0)){llik[i]<-log((1-pnorm(XB))+ pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
                else if((y[i]==1)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[3]-ZG)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(tau[2]-ZG)),mean=c(0,0),sigma=negsigma))}
                else if((y[i]==2)) {llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB,(ZG-tau[3])),mean=c(0,0),sigma=sigma))}     
		}
        llik<--1*sum(llik)					      
    	return(llik)
}	

############################################################################
#give the ZiOPC model resonable starting parameters
b<-rbind( -1.31, .32, 12.5, -.21, -.3, -1.01, -1.01, -1.01, 0.001603, -0.095935, 0.273732, 1.900761,-0.562974,-.4)

#optimize the ZiOPC model
output.ZIOPC<-nlm(f=ZIOPC, p=b, iterlim=500, data=BP.data, hessian=TRUE)

#show the output to get the coefficient estimates
print(output.ZIOPC)

#recover the variance covariance matrix
vcv<-solve(output.ZIOPC$hessian)

#show the variance covariance matrix to recover the standard errors
vcv 


############################################################################
#############################AIC Statistics#################################
############################################################################

#Ordered Probit
AIC.OP<--2*(-1432)+2*(7)
AIC.OP

#ZiOP
AIC.ZIOP<--2*(-1386)+2*(10)
AIC.ZIOP

#Small ZiOPC
AIC.ZIOPC1<--2*(-1374)+2*(11)
AIC.ZIOPC1

#Large ZiOPC
AIC.ZIOPC2<--2*(-1370)+2*(14)
AIC.ZIOPC2


############################################################################
###############################Vuong Tests##################################
############################################################################

#First, compare smaller ZiOP(C) models to OP model

#set the parameters
y<-BP.data[,8]			#political violence                                               
x1<-BP.data[,4]			#lgdppc                                           
x2<-BP.data[,9]			#parliament
x3<-BP.data[,7]			#disaster
x4<-BP.data[,5]			#oil
x5<-BP.data[,6]			#primary
x<-cbind(1,x1,x2)		#inflation stage covariates 		     	      
z<-cbind(x1,x2,x3,x4,x5)	#outcome stage covariates
n<-nrow(BP.data)		#number of observations
     
#Vuong test comparing small ZIOP to OP
m<-matrix(0,nrow=n,ncol=1)
  for(i in 1:n){					              					      
	ZG.OP<-(-0.2123*z[i,1]+-0.5380*z[i,2]+0.2203*z[i,3]+0.9073*z[i,4]+-0.4266*z[i,5])
	cut1.OP<--1.0729
	cut2.OP<--1.0729+exp(-0.1711)
	XB.ZIOP<-(18.78179*x[i,1]+-2.08193*x[i,2]+-0.29258*x[i,3])
	ZG.ZIOP<-(0.04125*z[i,1]+-0.09508*z[i,2]+0.26499*z[i,3]+1.70694*z[i,4]+-0.42220*z[i,5])
	cut1.ZIOP<-0.77186
	cut2.ZIOP<-0.77186+exp(-0.09820) 
        if((y[i]==0)){m[i]<-log(pnorm(cut1.OP-ZG.OP)/((1-pnorm(XB.ZIOP))+(pnorm(XB.ZIOP))*(pnorm(cut1.ZIOP-ZG.ZIOP))))}
        else if((y[i]==1)){m[i]<-log((pnorm(cut2.OP-ZG.OP)-pnorm(cut1.OP-ZG.OP))/((pnorm(XB.ZIOP))*((pnorm(cut2.ZIOP-ZG.ZIOP)) - (pnorm(cut1.ZIOP-ZG.ZIOP)))))}
        else if((y[i]==2)) {m[i]<-log((1-pnorm(cut2.OP-ZG.OP))/((pnorm(XB.ZIOP))*((1-pnorm(cut2.ZIOP-ZG.ZIOP)))))}     
	}
diff.m.sqrd<-(m-mean(m))^2
sum.d.m.s<-sum(diff.m.sqrd)
Vuong.OP.ZIOP<-((sqrt(n))*(1/n)*sum(m))/sqrt((1/n)*(sum.d.m.s))
Vuong.OP.ZIOP


#Vuong test comparing small ZIOPC to  OP
m<-matrix(0,nrow=n,ncol=1)
  for(i in 1:n){					              					      
	ZG.OP<-(-0.2123*z[i,1]+-0.5380*z[i,2]+0.2203*z[i,3]+0.9073*z[i,4]+-0.4266*z[i,5])
	cut1.OP<--1.0729
	cut2.OP<--1.0729+exp(-0.1711)
	XB.ZIOPC<-(11.6098*x[i,1]+-1.2809*x[i,2]+-0.3676*x[i,3])
	ZG.ZIOPC<-(0.3307*z[i,1]+0.3096*z[i,2]+0.1975*z[i,3]+1.1828*z[i,4]+-0.2367*z[i,5])
	cut1.ZIOPC<-2.7558
	cut2.ZIOPC<-2.7558+exp(-0.2139)
	rho<--0.8890          
	negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)
	sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
        if((y[i]==0)){m[i]<-log(pnorm(cut1.OP-ZG.OP)/((1-pnorm(XB.ZIOPC))+pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==1)){m[i]<-log((pnorm(cut2.OP-ZG.OP)-pnorm(cut1.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut2.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==2)) {m[i]<-log((1-pnorm(cut2.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(ZG.ZIOPC-cut2.ZIOPC)),mean=c(0,0),sigma=sigma)))}     
	}
diff.m.sqrd<-(m-mean(m))^2
sum.d.m.s<-sum(diff.m.sqrd)
Vuong.OP.ZIOPC<-((sqrt(n))*(1/n)*sum(m))/sqrt((1/n)*(sum.d.m.s))
Vuong.OP.ZIOPC



#Next, compare larger ZiOPC model to OP model

y<-BP.data[,8]			#political violence                                                    
x1<-BP.data[,4]			#lgdppc                                           
x2<-BP.data[,9]			#parliament
x3<-BP.data[,7]			#disaster
x4<-BP.data[,5]			#oil
x5<-BP.data[,6]			#primary
x<-cbind(1,x1,x2,x3,x4,x5)	#inflation stage covariates		     	      
z<-cbind(x1,x2,x3,x4,x5)	#outcome stage covariates
                 
#Vuong test comparing large ZIOPC to  OP
m<-matrix(0,nrow=n,ncol=1)
  for(i in 1:n){					              					      
	ZG.OP<-(-0.2123*z[i,1]+-0.5380*z[i,2]+0.2203*z[i,3]+0.9073*z[i,4]+-0.4266*z[i,5])
	cut1.OP<--1.0729
	cut2.OP<--1.0729+exp(-0.1711)
	XB.ZIOPC<-(12.43549*x[i,1]+-1.33011*x[i,2]+-0.02093*x[i,3]+-0.12895*x[i,4]+-1.91649*x[i,5]+0.85427*x[i,6])
	ZG.ZIOPC<-(0.20376*z[i,1]+-0.01334*z[i,2]+0.24377*z[i,3]+1.89114*z[i,4]+-0.57963*z[i,5])
	cut1.ZIOPC<-1.90998
	cut2.ZIOPC<-1.90998+exp(-0.21296)
	rho<--0.91232         
	negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)
	sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
        if((y[i]==0)){m[i]<-log(pnorm(cut1.OP-ZG.OP)/((1-pnorm(XB.ZIOPC))+pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==1)){m[i]<-log((pnorm(cut2.OP-ZG.OP)-pnorm(cut1.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut2.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(cut1.ZIOPC-ZG.ZIOPC)),mean=c(0,0),sigma=negsigma)))}
        else if((y[i]==2)) {m[i]<-log((1-pnorm(cut2.OP-ZG.OP))/(pmvnorm(lower=c(-Inf,-Inf), upper=c(XB.ZIOPC,(ZG.ZIOPC-cut2.ZIOPC)),mean=c(0,0),sigma=sigma)))}     
	}
diff.m.sqrd<-(m-mean(m))^2
sum.d.m.s<-sum(diff.m.sqrd)
Vuong.OP.ZIOPC<-((sqrt(n))*(1/n)*sum(m))/sqrt((1/n)*(sum.d.m.s))
Vuong.OP.ZIOPC

############################################################################
######################Likelihood Ratio Tests################################
############################################################################

#likelihood ratio test comparing the two ZiOPC models
l.ratio<-2*(-2768-(-2770))
pchisq(l.ratio, 3)





